/** 
 * @file test_gcomplex.cpp
 * @brief Complex number testing
 * @author Anshul Nigham
 * @date 2006-11-02
 */

#include <iostream>
#include <cmath>
#include <gsl/gsl_math.h>
#include "gcomplex.h"
#include "gcomplex_test.h"

using std::cout;

namespace gslcpp {

	/* TEMPLATE
   bool test_gcomplex_()
   {
   bool result = true;
   cout<<"testing gslcpp::gcomplex";
   gcomplex cplx;

   return result;
   }
	   */

	bool are_equal(gcomplex a, gsl_complex b)
	{
		return (FCMP(a.real(), GSL_REAL(b)) && FCMP(a.imag(), GSL_IMAG(b)));
	}

	bool test_gcomplex_defconstr()
	{
		bool result = true;
		gcomplex cplx;
		cout<<"testing gslcpp::gcomplex default constructor... ";
		result = result && cplx.real() == 0.0 && cplx.imag() == 0.0;
		return result;
	}

	bool test_gcomplex_rectconstr()
	{
		bool result = true;
		cout<<"testing gslcpp::gcomplex rect coord constructor... ";
		gcomplex cplx = gcomplex::rect(2.0,4.0);
		result = result && cplx.real() == 2.0 && cplx.imag() == 4.0;

		return result;
	}

	bool test_gcomplex_polarconstr()
	{
		bool result = true;
		cout<<"testing gslcpp::gcomplex polar coord constructor... ";
		double theta = M_PI/4.0;
		gcomplex cplx = gcomplex::polar(2.0,theta);
		result = result && FCMP(cplx.real(), 2.0*cos(theta)) &&
			FCMP(cplx.imag(), 2.0*sin(theta));

		return result;
	}

	bool test_gcomplex_set()
	{
		bool result = true;
		cout<<"testing gslcpp::gcomplex set function... ";
		gcomplex cplx;
		cplx.set(2.0,4.0);
		result = result && FCMP(cplx.real(), 2.0) && FCMP(cplx.imag(), 4.0);

		return result;
	}

	bool test_gcomplex_set_real()
	{
		bool result = true;
		cout<<"testing gslcpp::gcomplex set_real... ";
		gcomplex cplx;
		cplx.set_real(2.0);
		result = result && FCMP(cplx.real(), 2.0);

		return result;
	}

	bool test_gcomplex_set_imag()
	{
		bool result = true;
		cout<<"testing gslcpp::gcomplex set_imag... ";
		gcomplex cplx;
		cplx.set_imag(2.0);
		result = result && FCMP(cplx.imag(), 2.0);

		return result;
	}

	bool test_gcomplex_arg()
	{
		bool result = true;
		cout<<"testing gslcpp::gcomplex arg... ";
		gcomplex cplx;
		cplx.set(2.0,2.0);
		gsl_complex z = gsl_complex_rect(2.0,2.0);
		result = result && FCMP(gsl_complex_arg(z), cplx.arg());

		return result;
	}

	bool test_gcomplex_abs()
	{
		bool result = true;
		cout<<"testing gslcpp::gcomplex abs... ";
		gcomplex cplx;
		cplx.set(2.0,2.0);
		gsl_complex z = gsl_complex_rect(2.0,2.0);
		result = result && FCMP(gsl_complex_abs(z), cplx.abs());

		return result;
	}

	bool test_gcomplex_abs2()
	{
		bool result = true;
		cout<<"testing gslcpp::gcomplex abs2... ";
		gcomplex cplx;
		cplx.set(2.0,2.0);
		gsl_complex z = gsl_complex_rect(2.0,2.0);
		result = result && FCMP(gsl_complex_abs2(z), cplx.abs2());

		return result;
	}

	bool test_gcomplex_logabs()
	{
		bool result = true;
		cout<<"testing gslcpp::gcomplex logabs... ";
		gcomplex cplx;
		cplx.set(2.0,2.0);
		gsl_complex z = gsl_complex_rect(2.0,2.0);
		result = result && FCMP(gsl_complex_logabs(z), cplx.logabs());

		return result;
	}

	bool test_gcomplex_operators()
	{
		bool result = true;
		cout<<"testing gslcpp::gcomplex operators... ";
		gcomplex c1, c2, cplus, cminus, cmul, cdiv;
		c1.set(2.0,5.0); c2.set(5.0,2.0);

		gsl_complex z1 = gsl_complex_rect(2.0,5.0);
		gsl_complex z2 = gsl_complex_rect(5.0,2.0);
		gsl_complex zplus,zminus,zmul,zdiv;

		// test complex (op) complex

		cplus = c1+c2; zplus = gsl_complex_add(z1,z2);
		cminus = c1-c2; zminus = gsl_complex_sub(z1,z2);
		cmul = c1*c2; zmul = gsl_complex_mul(z1,z2);
		cdiv = c1/c2; zdiv = gsl_complex_div(z1,z2);

		result = result && are_equal(cplus,zplus) &&
			are_equal(cminus, zminus) &&
			are_equal(cmul, zmul) &&
			are_equal(cdiv, zdiv);

		// test complex (op) real
		double r = 5.0;

		cplus = c1+r; zplus = gsl_complex_add_real(z1,r);
		cminus = c1-r; zminus = gsl_complex_sub_real(z1,r);
		cmul = c1*r; zmul = gsl_complex_mul_real(z1,r);
		cdiv = c1/r; zdiv = gsl_complex_div_real(z1,r);

		result = result && are_equal(cplus,zplus) &&
			are_equal(cminus, zminus) &&
			are_equal(cmul, zmul) &&
			are_equal(cdiv, zdiv);

		cplus = c1.add_imag(r); zplus = gsl_complex_add_imag(z1,r);
		cminus = c1.sub_imag(r); zminus = gsl_complex_sub_imag(z1,r);
		cmul = c1.mul_imag(r); zmul = gsl_complex_mul_imag(z1,r);
		cdiv = c1.div_imag(r); zdiv = gsl_complex_div_imag(z1,r);

		result = result && are_equal(cplus,zplus) &&
			are_equal(cminus, zminus) &&
			are_equal(cmul, zmul) &&
			are_equal(cdiv, zdiv);

		return result;
	}

	bool test_gcomplex_trig()
	{
		gcomplex c;
		c.set(1.5,2.5);
		gsl_complex z = gsl_complex_rect(1.5,2.5);
		bool result = true;
		
		cout<<"testing gslcpp::gcomplex trigonometric functions... ";

		result = result && are_equal(c.sin(), gsl_complex_sin(z)) &&
			are_equal(c.cos(), gsl_complex_cos(z)) &&
			are_equal(c.tan(), gsl_complex_tan(z)) &&
			are_equal(c.sec(), gsl_complex_sec(z)) &&
			are_equal(c.csc(), gsl_complex_csc(z)) &&
			are_equal(c.cot(), gsl_complex_cot(z));

		return result;
	}

	bool test_gcomplex_invtrig()
	{
		gcomplex c;
		c.set(1.5,2.5);
		gsl_complex z = gsl_complex_rect(1.5,2.5);
		double r = 2.2;
		bool result = true;
		
		cout<<"testing gslcpp::gcomplex inverse trigonometric functions... ";

		result = result && are_equal(c.arcsin(), gsl_complex_arcsin(z)) &&
			are_equal(c.arccos(), gsl_complex_arccos(z)) &&
			are_equal(c.arctan(), gsl_complex_arctan(z)) &&
			are_equal(c.arcsec(), gsl_complex_arcsec(z)) &&
			are_equal(c.arccsc(), gsl_complex_arccsc(z)) &&
			are_equal(c.arccot(), gsl_complex_arccot(z));

		result = result &&
			are_equal(gcomplex::arcsin(r), gsl_complex_arcsin_real(r)) &&
			are_equal(gcomplex::arccos(r), gsl_complex_arccos_real(r)) &&
			are_equal(gcomplex::arcsec(r), gsl_complex_arcsec_real(r)) &&
			are_equal(gcomplex::arccsc(r), gsl_complex_arccsc_real(r));
		
		return result;
	}

	bool test_gcomplex_hyper()
	{
		gcomplex c;
		c.set(1.5,2.5);
		gsl_complex z = gsl_complex_rect(1.5,2.5);
		bool result = true;
		
		cout<<"testing gslcpp::gcomplex hyperbolic functions... ";

		result = result && are_equal(c.sinh(), gsl_complex_sinh(z)) &&
			are_equal(c.cosh(), gsl_complex_cosh(z)) &&
			are_equal(c.tanh(), gsl_complex_tanh(z)) &&
			are_equal(c.sech(), gsl_complex_sech(z)) &&
			are_equal(c.csch(), gsl_complex_csch(z)) &&
			are_equal(c.coth(), gsl_complex_coth(z));

		return result;
	}

	bool test_gcomplex_invhyper()
	{
		gcomplex c;
		c.set(1.5,2.5);
		gsl_complex z = gsl_complex_rect(1.5,2.5);
		double r = 2.2;
		bool result = true;
		
		cout<<"testing gslcpp::gcomplex inverse hyperbolic functions... ";

		result = result && are_equal(c.arcsinh(), gsl_complex_arcsinh(z)) &&
			are_equal(c.arccosh(), gsl_complex_arccosh(z)) &&
			are_equal(c.arctanh(), gsl_complex_arctanh(z)) &&
			are_equal(c.arcsech(), gsl_complex_arcsech(z)) &&
			are_equal(c.arccsch(), gsl_complex_arccsch(z)) &&
			are_equal(c.arccoth(), gsl_complex_arccoth(z));

		result = result &&
			are_equal(gcomplex::arccosh(r), gsl_complex_arccosh_real(r)) &&
			are_equal(gcomplex::arctanh(r), gsl_complex_arctanh_real(r));
		
		return result;
	}

	bool gcomplex_runtest()
	{
		bool result = true;
		tfarray test_functions;
		test_functions.push_back(test_gcomplex_defconstr);
		test_functions.push_back(test_gcomplex_rectconstr);
		test_functions.push_back(test_gcomplex_polarconstr);
		test_functions.push_back(test_gcomplex_set);
		test_functions.push_back(test_gcomplex_set_real);
		test_functions.push_back(test_gcomplex_set_imag);
		test_functions.push_back(test_gcomplex_arg);
		test_functions.push_back(test_gcomplex_abs);
		test_functions.push_back(test_gcomplex_abs2);
		test_functions.push_back(test_gcomplex_logabs);
		test_functions.push_back(test_gcomplex_operators);
		test_functions.push_back(test_gcomplex_trig);
		test_functions.push_back(test_gcomplex_invtrig);
		test_functions.push_back(test_gcomplex_hyper);
		test_functions.push_back(test_gcomplex_invhyper);

		for (int i=0; i<test_functions.size(); ++i) {
			bool tmp_result = test_functions[i]();
			if (tmp_result) cout<<"PASS\n";
			else cout<<"FAIL!\n";
			result = result && tmp_result;
		}
		return result;
	}

}

